
#####
# Leader of the Pack?
# Changes in ‘Wolf Warrior Diplomacy’ After 
# A Politburo Collective Study Session

# Authors:
# Samuel Brazys
# Alexander Dukalskis
# Stefan Müller


# 03_plot_results.R creates the graphs for the paper and appendix
# based on the regression models run in Stata and LSS data
####

library(tidyverse) # CRAN v1.3.1
library(cowplot)   # CRAN v1.1.1
library(forcats)   # CRAN v0.5.1
library(here)      # CRAN v1.0.1

# import custom ggplot2 theme

theme_baser <- function (){
    theme_minimal()  %+replace%
        theme(panel.grid.minor.x = element_blank(),
              panel.grid.minor.y = element_blank(),
              panel.grid.major.x = element_blank(),
              panel.grid.major.y = element_blank(),
              panel.border = element_rect(fill = NA,colour = "black", size = 0.5,
                                          linetype = "solid"),
              legend.title = element_text(size = 15),
              plot.title = element_text(size = 15, face = "bold",
                                        vjust = 1.5, hjust = 0.5,
                                        margin=margin(5, 5, 5 ,5)),
              plot.caption = element_text(colour = "grey30", size = 11, hjust = 1),
              legend.position = "bottom",
              axis.ticks.y = element_line(size = 0.3),
              axis.ticks.x = element_line(size = 0.3),
              axis.ticks.length = unit(0.2, "cm"),
              legend.text=element_text(size = 13),
              panel.background = element_rect(fill='white'), #transparent panel bg
              plot.background = element_rect(fill='white', color= "white"), #transparent plot bg
              strip.text = element_text(size = 15, hjust = 0.5,
                                        face = "bold",
                                        margin = margin(b = 5, r = 5, l = 5, t = 5)),
              axis.text.y = element_text(colour = "black", size = 13,
                                         hjust = 1),
              axis.text.x = element_text(colour = "black", size = 13),
              axis.title = element_text(size = 13, hjust = 0.5))
}


# set theme
theme_set(theme_baser())


# load data with LSS scores
dat <- read_csv("data_tweets_en.csv")

table(dat$country)

dat_day <- dat %>% 
    mutate(pre_post = ifelse(date < as.Date("2021-06-01"), "Before Announcement", "After Announcement")) %>% 
    group_by(date, pre_post, type) %>% 
    summarise(mean_wolfy = mean(wolfy_score, na.rm = TRUE))

dat_day$pre_post <- relevel(factor(dat_day$pre_post), ref = "Before Announcement")


dat_day_all <- dat %>% 
    filter(type != "Media Account") %>% # exclude media accounts
    mutate(pre_post = ifelse(date < as.Date("2021-06-01"), "Before Announcement", "After Announcement")) %>% 
    group_by(date, pre_post) %>% 
    summarise(mean_wolfy = mean(wolfy_score, na.rm = TRUE),
              sd_wolfy = sd(wolfy_score),
              n = n()) %>%
    mutate(se_wolfy = sd_wolfy / sqrt(n),
           ci_lower_90 = mean_wolfy - 1.645 * se_wolfy,
           ci_upper_90 = mean_wolfy + 1.645 * se_wolfy,
           ci_lower_95 = mean_wolfy - 1.96 * se_wolfy,
           ci_upper_95 = mean_wolfy + 1.96 * se_wolfy)


ggplot(
    filter(dat_day_all, between(date, as.Date("2021-03-01"), as.Date("2021-09-01"))),
    aes(date, y = mean_wolfy)) +
    geom_linerange(aes(ymin = ci_lower_90, max = ci_upper_90), 
                   size = 1.05,
                   colour = "grey60") +
    geom_linerange(aes(ymin = ci_lower_95, max = ci_upper_95), 
                   colour = "grey60") +
    geom_point(size = 2, shape = 16, colour = "grey60")  +
    geom_smooth(colour = "black", method = "gam", n = 1000, se = FALSE) +
    geom_vline(xintercept = as.Date("2021-06-01"), 
               linetype = "dashed") +
    scale_fill_manual(values = c("grey70", "black")) +
    scale_x_date(date_labels = "%b %Y", date_breaks = "1 month") +
    labs(x = NULL, y = "Average WWD Score") +
    annotate("label", x = as.Date("2021-03-05"),
             y = 0.25, label = "Friendly",
             label.size = NA,
             colour = "grey30", size = 4) +
    annotate("text", x = as.Date("2021-05-26"),
             hjust = 1,
             y = -.5, label = "Before announcement",
             colour = "grey30", size = 4) +
    annotate("text", x = as.Date("2021-06-10"),
             hjust = 0,
             y = -.5, label = "After announcement",
             colour = "grey30", size = 4) +
    annotate("segment",
             x = as.Date("2021-03-15"),
             xend = as.Date("2021-03-15"),
             y = 0.2, yend = 0.3, colour = "grey30",
             size = 0.5,
             arrow = arrow(angle = 25, length = unit(0.3, "cm"))) +
    annotate("segment",
             x = as.Date("2021-05-26"),
             xend = as.Date("2021-04-15"),
             y = -0.55, yend = -0.55, colour = "grey30",
             size = 0.5,
             arrow = arrow(angle = 25, length = unit(0.3, "cm"))) +
    annotate("segment",
             x = as.Date("2021-06-10"),
             xend = as.Date("2021-07-21"),
             y = -0.55, yend = -0.55, colour = "grey30",
             size = 0.5,
             arrow = arrow(angle = 25, length = unit(0.3, "cm"))) +
    theme(legend.position = "none")
ggsave("fig_01.pdf",
       width = 9, height = 5)

# change scipen option to add nicer labels to plots
options(scipen=200000000)

# function to plot the diff-in-diff estimates for varying windows
plot_window <- function(x) {
    
    # make sure these are the correct coefficients
    tab_windows_clean <- x[6:7, ]
    
    # diff needs to be included
    table(tab_windows_clean$X)
    
    print(head(tab_windows_clean$X))
    
    dat_windows_long_coef <- tab_windows_clean[1, ] %>%
        select(-X) %>% 
        gather(days, coef) %>% 
        mutate(coef = as.numeric(coef))
    
    dat_windows_long_se <- tab_windows_clean[2, ] %>% 
        select(-X) %>% 
        gather(days, se) %>% 
        mutate(se = as.numeric(se))
    
    dat_windows_long <- left_join(dat_windows_long_coef, 
                                  dat_windows_long_se,
                                  by = "days") %>% 
        mutate(days = str_replace_all(days, "all_", " ")) %>% 
        mutate(days = as.numeric(days))
    
    ggplot(dat_windows_long, 
           aes(x = days, y = coef)) +
        geom_linerange(aes(ymin = coef - 1.645 * se,
                           ymax = coef + 1.645 * se),
                       size = 1) +
        geom_linerange(aes(ymin = coef - 1.96 * se,
                           ymax = coef + 1.96 * se),
                   size = 0.5) +
    geom_point(size = 1.5) +
        geom_hline(yintercept = 0, size = 0.8,
                   linetype = "dashed", colour = "red") +
    scale_x_continuous(breaks = c(seq(1, 130, 10))) +
    annotate("label", x = 16,
             hjust = 0,
             y = 0.45, 
             label = "OECD accounts friendlier",
             label.size = NA,
             colour = "grey30", size = 4.5) +
    scale_y_continuous(limits = c(-0.6, 0.7),
                       breaks = c(-0.5, -0.25, 0,  0.25, 0.5)) +
    annotate("segment",
             x = 14,
             xend = 14,
             y = 0.35, yend = 0.55, colour = "grey30",
             size = 0.5,
             arrow = arrow(angle = 25, length = unit(0.3, "cm"))) +
    labs(x = "Window of Days Around Announcement", y = "Difference OECD vs Non-OECD") +
        theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 17,
                                        margin=margin(b=5, r = 5, l = 5, t = 5)))
}




# function for main plot
plot_window_main <- function(x) {
    
    # make sure these are the correct coefficients
    tab_windows_clean <- x[6:7, ]
    
    # diff needs to be included
    table(tab_windows_clean$X)
    
    print(head(tab_windows_clean$X))
    
    dat_windows_long_coef <- tab_windows_clean[1, ] %>%
        select(-X) %>% 
        gather(days, coef) %>% 
        mutate(coef = as.numeric(coef))
    
    dat_windows_long_se <- tab_windows_clean[2, ] %>% 
        select(-X) %>% 
        gather(days, se) %>% 
        mutate(se = as.numeric(se))
    
    dat_windows_long <- left_join(dat_windows_long_coef, 
                                  dat_windows_long_se,
                                  by = "days") %>% 
        mutate(days = str_replace_all(days, "all_", " ")) %>% 
        mutate(days = as.numeric(days))
    
    ggplot(dat_windows_long, 
           aes(x = days, y = coef)) +
        geom_linerange(aes(ymin = coef - 1.645 * se,
                           ymax = coef + 1.645 * se),
                       size = 1) +
        geom_linerange(aes(ymin = coef - 1.96 * se,
                           ymax = coef + 1.96 * se),
                       size = 0.5) +
        geom_point(size = 1.5) +
        geom_hline(yintercept = 0, 
                   linetype = "dashed", colour = "red",
                   size = 0.8) +
        scale_x_continuous(breaks = c(seq(1, 130, 10))) +
        annotate("label", x = 16,
                 hjust = 0,
                 y = 0.3, 
                 label = "OECD accounts friendlier",
                 label.size = NA,
                 colour = "grey30", size = 4.5) +
        scale_y_continuous(limits = c(-0.45, 0.45),
                           breaks = c(-0.4, -0.2, 0, 0.2, 0.4)) +
    annotate("segment",
             x = 15,
             xend = 15,
             y = 0.2, yend = 0.4,
             colour = "grey30",
             size = 0.4,
             arrow = arrow(angle = 25, length = unit(0.3, "cm"))) +
        labs(x = "Window of Days Around Announcement", y = "Difference OECD vs Non-OECD") 

}


# import results from Stata do file and create plots

tab_windows <- read.csv("window_results.csv")


plot_window_main(tab_windows)
ggsave("fig_02.pdf",
       width = 9, height = 5)

# institutions
tab_windows_inst <- read.csv("window_results_INST.csv")

p_inst <- plot_window(tab_windows_inst) +
    labs(title = "Institutions")

# ambassadors
tab_windows_amb <- read.csv("window_results_AMB.csv")

p_amb <- plot_window(tab_windows_amb) +
    labs(title = "Ambassadors")

# staff
tab_windows_stf <- read.csv("window_results_STF.csv")

p_stf <- plot_window(tab_windows_stf) +
    labs(title = "Staff")


# diplomats
tab_windows_dip <- read.csv("window_results_DIP.csv")

p_dip <- plot_window(tab_windows_dip) +
    labs(title = "Diplomats")

# combine plots using the plot_grid function from cowplot
cowplot::plot_grid(p_amb, p_dip, p_inst, p_stf, nrow = 2)
ggsave("fig_a01.pdf",
       width = 14, height = 10)



# pre-trend plots

dat_pretrend_all <- read.csv("pretrend_all.csv",
                             na.strings = "") %>% 
    mutate(treat = ifelse(treat == "Treated", "OECD Accounts\n(Treated)",
                          ifelse(treat == "Not Treated", "Non-OECD Accounts\n(Not Treated)", NA)))

dat_pretrend_all$treat <- forcats::fct_rev(factor(dat_pretrend_all$treat))

dat_pretrend_all$accounts <- "All (Diplomats and Media)"


dat_pretrend_amb <- read.csv("pretrend_amb.csv",
                             na.strings = "") %>% 
    filter(!is.na(treat)) %>% 
    mutate(treat = ifelse(treat == "Treated", "OECD Accounts\n(Treated)",
                          ifelse(treat == "Not Treated", "Non-OECD Accounts\n(Not Treated)", NA)))

dat_pretrend_amb$treat <- forcats::fct_rev(factor(dat_pretrend_amb$treat))

dat_pretrend_amb$accounts <- "Ambassadors"


dat_pretrend_dip <- read.csv("pretrend_dip.csv",
                             na.strings = "") %>% 
    filter(!is.na(treat)) %>% 
    mutate(treat = ifelse(treat == "Treated", "OECD Accounts\n(Treated)",
                          ifelse(treat == "Not Treated", "Non-OECD Accounts\n(Not Treated)", NA)))



dat_pretrend_dip$treat <- forcats::fct_rev(factor(dat_pretrend_dip$treat))

dat_pretrend_dip$accounts <- "Diplomats"



dat_pretrend_int <- read.csv("pretrend_int.csv",
                             na.strings = "") %>% 
    filter(!is.na(treat)) %>% 
    mutate(treat = ifelse(treat == "Treated", "OECD Accounts\n(Treated)",
                          ifelse(treat == "Not Treated", "Non-OECD Accounts\n(Not Treated)", NA)))


dat_pretrend_int$treat <- forcats::fct_rev(factor(dat_pretrend_int$treat))

dat_pretrend_int$accounts <- "Institutions"


dat_pretrend_stf <- read.csv("pretrend_stf.csv",
                             na.strings = "") %>% 
    filter(!is.na(treat)) %>% 
    mutate(treat = ifelse(treat == "Treated", "OECD Accounts\n(Treated)",
                          ifelse(treat == "Not Treated", "Non-OECD Accounts\n(Not Treated)", NA)))


dat_pretrend_stf$treat <- forcats::fct_rev(factor(dat_pretrend_stf$treat))

dat_pretrend_stf$accounts <- "Staff"



dat_pretrend_all$accounts <- "All"


dat_pretends <- bind_rows(dat_pretrend_all,
                          dat_pretrend_amb,
                          dat_pretrend_dip,
                          dat_pretrend_int,
                          dat_pretrend_stf)


ggplot(dat_pretends, aes(x = as.numeric(day_c),
                         y = WDD, colour = time)) +
    geom_vline(xintercept = 0, linetype = "dashed", colour = "red") +
    geom_smooth(method = "loess") +
    scale_colour_manual(values = c("black", "grey50")) +
    facet_grid(accounts~treat) +
    scale_x_continuous(breaks = c(seq(-25, 25, 5))) +
    coord_cartesian(ylim = c(-1, 1)) +
    labs(x = "Days around Announcement", y = "WDD Score") +
    theme(legend.position = "none", strip.text.y = element_text(angle = 0, hjust = 0)) 
ggsave("fig_a02.pdf", 
       width = 9, height = 6)
